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INTRODUCTION 


Thu orbits of the satellites of the outer 
planets are poorly known, due to lack of 
attention over the past half century. We 
have been developing a new theory of Saturn's 
satellites Enceladus and Dione which is 
literal (all constants of integration 
appear explicitly), canonically invariant 
(the Hori-Lie method is used), and which 
correctly handles the eccentricity- type 
resonance between the two satellites. The 
algebraic manipulations are being performed 
using the TRIGMAN formula manipulation 
language, and the programs have been 
developed so that with minor modifications 
they can be used on the Mimas-Tethys 
and Titan-Hyper ion systems. 



Our current knowledge about the orbits of the outer planets' satellites 

is quite poor. The theories now in use date back to the early part of this 

century, and no longer provide accurate ephemerides. At the same time, 

these satellites have been virtually unobserved for the better part of this 

1 

century. As a result, it has been estimated by Mulholland that the errors 
in the predicted coordinates may range into the tens of thousands of kilometers. 


A primary purpose of the upcoming missions to the outer planets will be to 
determine Improved positions for these satellites, in addition to making 
physical observations. It Is therefore very important that their orbits be 
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rodisoussed in the light of modern as well as older observations. The re- 
observation of the satellites is already well under way at McDonald 
Observatory , and to go along with it, a new theoretical discussion of 
several of the satellites of Saturn has been undertaken. The system 
Enceladus-Dione is the first for which results have been obtained, and 
it is on this system that we report in this paper, 

REMARKS ON PERTURBATION THEORY 

Traditionally, there have been two principal approaches to solving perturbed 
motion in Dynamical Astronomy. One, numerical integration of the equations 
of motion, has much to commend it, and in fact, is the usual approach these 
days. It is straightforward and easy to implement, and can model a wide 
variety of situations accurately. The other approach, the so-called 
"literal" method, attempts to produce explicit formulae for the determina- 
tion of a body's position at any arbitrary time. It is difficult to 
implement, as it requires special series manipulation programs — programs 
which are very complicated to write and which require unfamiliar techniques 
in Iheir construction. In addition, numerous technical problems arise in 
practice, the most common being the ease with which it is possible to 
build up very long expressions which threaten to overrun the entire 
core memory available. 

Despite these disadvantages of the literal method, there are several 
important advantages. One is that the formulae may be good for very long 
periods of time, and it takes no more effort to obtain a prediction for a 
time far removed from the epoch of the theory than it takes to determine a 
prediction for a time near the epoch. With the numerical approach this 
is not so. A second advantage is that the constants of integration, and 
other parameters of the theory (e.g. , the masses and oblateness parameters) 
can be made to appear explicitly in the formulae. This means that the 
theory can be used even if better values of the parameters are determined. 
Yet a third advantage is that each particular term in a formula arises in 
a predictable way and has a predictable effect on the motion of the bodies. 
This gives one deeper insight into the problem, particularly when resonances 
are of importance fas they are in this case). Thus, there are good reasons 
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for the interest which has been shown by celestial mechanicians in the 
recent progress which has been made in the area of computerized algebra 
'or dynamical astronomy. 

SYSTEMS FOR COMPUTERIZED ALGEBRA 

In order to see what is presently possible, let us look a little more 

deeply into the kinds of processors which are currently available. We 

will concentrate on the Poisson series processor, which is the most commonly 
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available. Such programs have been written by Barton , Broucke , Brumberg , 
6 7 8 

Chapront et al, , Jefferys and Rom . With several variations on a theme, 
they are capable of adding, subtracting, multiplying, differentiating, 
integrating, and performing certain substitution operations on series 
of the general form of Eq, (1) : 

S = E •£ C k (x) cos ( k. y ) + S k (x) sin ( k* y ) } , (1) 

where k = (k^, ...,k n ) is a vector of integers (positive, negative and 
zero), x = (x^, ...,x m ) and y = (y^, ...y n ) are vectors of polynomial 
variables and trigonometric variables, respectively, C^(x) and S^(x) are 
polynomials in the x variables having rational or real coefficients 
(depending on the system), and k»y = k^y 1 + ... + l< n y n is the dot 
product. Thus a Poisson series is a multiple Fourier series in a linear 
combination (with integer multipliers) of trigonometric variables, 
with its coefficients drawn from a ring of polynomials in several 
independent variables. 

Series of the form of Eg. (1) are of particular importance in dynamical 
astronomy because the motion of bodies in the solar system is in general 
quasiperiodic, and can be expressed in the form of Eq. (1) provided that 
the angles are linear functions of the time (with, in general, frequencies 
that are not rational ratios of each other) . More generally, one can 
write the equations of motion in the form 

x = e f(x,y) ^ 

y = v + sg(x,y) 

where the components of the vector functions f and g are Poisson series, 

V Is the vector of the mean motions of the angles, and 0 is a small 
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quantity. To first order we can write down the formal solution of these 
equations as 


x = x Q + e Jf (x 0 ,y Q -Put)dt, 
y = y Q + Vt + e [ g (x 0 ,y 0 +vt)dt, 


( 3 ) 


calculations which can easily be performed with any of the Poisson series 
processors currently available. 


RESONANCES 


A difficulty arises - one which is of great importance in the case of 
Enceladus and Dione - whenever there exi.sts a nearly exact relationship 
between the components of the vector v such that, for some vector of 
integers j, one has v.j w 0. Since, in general, all possible vectors of 
integers j are in fact represented in the series f and g, and since any 
real number can be approximated arbitrarily closely by a ratio of integers, 
this situation (as has been known for a long time) seriously threatens the 
convergence of Eqs. (3) and their higher order approximations. What 
happens is that in taking the integrals in Eqs. (3), divisors of the form 
v*;j appear. Since these divisors are dense on the real line, and 
in particular possess 0 as a limit point, the convergence of the series 
is in question. 

It is usual to ignore this problem as being of little practical interest 
(although of great theoretical interest) , unless one of the terms having 
such a "small divisor" also happens to have a large coefficient. This is 
indeed the case for Enceladus and Dione, and is due to the fact that the 
mean motion of Enceladus is almost exactly twice that of Dione. In fact, 
there are two terms which are greatly amplified by the presence of a 
small divisor: one is the term in the angular argument X “ ~ 2\Jj ' +'D3 / , 

where is the longitude of Enceladus, *C5 is the longitude of its perisaturnium, 
and the primed variables represent the same quantities for Dione (Fig. 1). 

This argument has a mean motion of 30° per yeai' , while most terms have mean 
motions which are over 1000 times as great. As a result, this term is 
amplified in Eqs. (3) by a factor which is over 1000 times greater than 
that for a typical term. In addition, the unamplified coefficient of this 
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term is one of the largest in the entire series. 

An even more serious problem is posed by the term in the argument 
c = (!j - 2i!|' + TU , for which the mean motion averages out to precisely zero, 
a situation which Garfinkel^ has called "deep resonance." The presence 
of an exact resonance such as this one causes the character of the solution 
to change completely. The angle 8, instead of either increasing or decreasing 
at a more or less constant, albeit slow rate, librates , or oscillates 
about some mean value (here, 0°). The reason for this behavior is that 
the phase plane turns out, in the lowest order, to be the same as that 
for a nonlinear pendulum (Fig. 2). Because of this, a different approach 
to this problem must be taken than the naive one indicated in Eqs. (3). 

DESCRIPTION OF THE PRESENT THEORY 

The detailed equations employed in the present theory will be published 
elsewhere"^, and therefore will not be given here. Instead, we shall try to 
give an overall view of the theory and point out some of its more interest- 
ing Features. 

Our approach has been to make use of one of the recently developed Lie 
series methods for canonical perturbation theory (since the present problem 
Is conservative and can be derived from a Hamiltonian). The Hori-Lie 
method that we have used has a number of advantages from our point of view: 

(1) It is canonically invariant, so that the theory, once completed, 
is valid for any coordinate system. The reason for this is 
that the Poisson bracket, on which the method is based, is 

a coordinate- independent notion. 

(2) It is readily adapted for systematic calculation by machine. 

(3) The perturbations in any desired quantity, including directly 
observed quantities such as the position or velocity vector, 
are obtained by application of a simple and uniform formula 
which is readily programmed, 

(tr) Perturbations are derived from the reduced Hamiltonian 

and only one other single function, the determining function. 

Since the actual coordinate system chosen is of little importance in applying 
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the method, it can be chosen for convenience. We have chosen to work in, 
variables closely related to the classical elements of semimajor axis, 
eccentricity and sine inclination, since it is relatively straightforward 
to write down the defining formulas for the disturbing function in terms 
of them. One modification we have adopted is to use as our actual 
variables, deviations from nominal values; the nominal values have been 
adopted from Struve’s theory. The advantage of this approach is that only 
low order polynomials in the deviations need to be retained in the develop- 
ment of the Hamiltonian in order to retain sufficient accuracy. Thus the 
theory can be made more compact. 

The first, and most difficult task is to develop the term in the 
Hamiltonian that describes the interaction between the two satellites. 
Because the ratio of the semimajor axes of the two satellites is so large 
(about 0.63), the series for this term converges very slowly. It is 
relatively straightforward to obtain an expression for the square of the 
distance between the two bodies, and a very short series results. However, 
the Newtonian interaction is proportional to the inverse of the distance, 
and is developed in powers of the ratio 0.63, which only slowly approaches 
zero. 

Many methods have been developed for performing the calculation, most of 

them prior to the advent of modern Poisson series processors. Our method, 

which harks back to one described by Brouwer and Clemen ne , has been 

simply to raise the expression for the square of the mutual distance to the 

inverse one-half power. We have employed for this an iterative technique 

14 

advocated by Broucke which converges quadratically, i. e., the error at 
each step is of the order of the square of the error at the preceding step. 
This method has been very satisfactory and converges with great rapidity. 
The rest of the Hamiltonian is very simple to calculate, having none of 
the problems of convergence associated with the above term. This includes 
the "indirect” term in the Newtonian potential which arises from choosing 
the origin of the system at the center of Saturn, and also the terms which 
arise from Saturn’s gravitational quadrupole moment. 

Our next step is to perform a coordinate transformation, using the Hori-Lie 
method, to eliminate from the Hamiltonian all terms in angular combinations 
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which arc not simple multiples of the resonant argument 0. This includes 
the terms in the shallowly resonant argument which, although they are 
considerably amplified and contribute significantly to the solution, 
nevertheless do not change the actual character of the solution as does 
the argument 0. A determining function is obtained which allows the 
calculation of the effects of all the angular arguments except 0 on the 
motion of the two satellites. We have also programmed a subroutine which 
allows the perturbations in any desired quantity to be determined, insofar 
as they are affected by these angular arguments. 

The remaining system, which contains the deeply resonant argument 0, must 
now be solved. Since this is the only angular argument remaining, the 
system has effectively been reduced to one degree of freedom. Looking 
again at Fig. 2, it can be seen that the motion in 0 and e is basically 
oscillatory in nature. In the actual case, the amplitude of the motion 
is quite small, only about 1H T of arc in 0, so that it is not even necessary 
to worry about the nonlinear nature of the pendulum presented when the 
amplitude is large. In the latter case, elliptic integrals are required 
to get the solution; but in our case, it is sufficient to approximate 
the system with a linear harmonic oscillator. This we have done by intro- 
ducing an appropriate new set of canonical elements, centered on the 
origin in Fig. 2, such that in these coordinates the equations of motion 
are reduced to an (amplitude dependent:) harmonic oscillator. The solution 
to this being well-known, we have therefore arrived at the complete solution 
of the problem. 

FINAL FORM OF THE THEORY 

In this form, the theory is not immediately useful, particularly if one is 
lacking a Poisson series processor. Therefore we have employed another method 
in order to make the theory accessible to the "outside world." Since 
FORTRAN is nearly ubiquitous as a programming language, we have chosen to 
have our Poisson series processor punch out FORTRAN- compatible expressions 
for each of the formulae required to obtain perturbations in any desired 
quantity. Thus, given numerical values of the various parameters in the 


theory, these expressions can be used to compute numerical values of any 
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desired quantity at any time. 

We are in the process of setting up a series of subroutines for the computation 
of all the necessary quantities. As can be imagined, many of the expressions 
are quite long and involved, so a fair amount of machine time is consumed 
in their evaluation. However, modern machines are really quite fast, 
and experience has shown that the time required for a single evaluation 
is not prohibitive (amounting to many evaluations per second). Thus, 
our method is quite practical. Furthermore, anyone who has access to a 
competent FORTRAN compiler will be able to use our theory. 

At the present time, we are setting up programs for the comparison of our 

literal theory with an ephemeris generated by numerical Integration, as 

a check on the accuracy of the theory. At the same timq, we are programming 

a set of subroutines to obtain definitive values for the parameters of 

the theory (constants of integration, masses and oblateness parameters) 

by a differential correction pi>ocedure. Since the data are still being 

gathered for this step, we are checking this phase also against numerical 

integrations, but eventually we expect to use both the old and new observations 

to produce a definitive theory. We will be greatly helped in this step 

15 

by the recent publication by Pierce of a bibliography of observational 
material on the satellites of Saturn. 

Eventually it is hoped to include observations of the satellites from the 
space missions themselves. To this end, we are working for an ultimate 
predictive capability for the theory of ±1 km. 

Finally, we mention that the programs we have written are adaptable to more 
problems than the Enceladus-Dione case. Both the Mimes-Tethys and Titan- 
Hyperion systems also have strong resonances (involving different angular 
arguments) . We have attempted to write the programs so that with minor 
changes they can be applied to these systems as well. 
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FIGURE CAPTIONS 

Figure 1: Geometry of the Orbits 

Figure 2: Resonance Phase Plane 



satellite 


perisaturnium 

orbit plane 
reference plane 


gj = >0, + cu 
v// = £1 + co + y 



